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Avalanche of Bifurcations and Hysteresis in a Model of Cellular Differentiation 
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Cellular differentiation in a developping organism is studied via a discrete bistable reaction- 
diffusion model. A system of undifferentiated cells is allowed to receive an inductive signal emenating 
from its environment. Depending on the form of the nonlinear reaction kinetics, this signal can 
trigger a series of bifurcations in the system. Differentiation starts at the surface where the signal is 
received, and cells change type up to a given distance, or under other conditions, the differentiation 
process propagates through the whole domain. When the signal diminishes hysteresis is observed. 

PACS numbers: 87.16.Ac, 87.17.Ee,87.18.Hf 



I. INTRODUCTION 

An adult higher organism, hke a human, has some hun- 
dreds of functionaUy different cell types. The genetic 
code stored by the DNA in the cell nucleus is identical 
in these cells. This potential information, however, is 
not utilized completely by the cells as many genes stay 
in a dormant, unexpressed state. The spectrum of genes 
which are expressed and functioning varies from cell type 
to cell type. One of the most fascinating questions in 
modern biology is how a certain cell or a group of cells 
finds its place and special task (cell type) in a developing 
organism. 1] 

From the point of view of dynamical systems, different 
cell types in the organism can be associated with different 
attractors of the common nonlinear internal dynamics of 
the cells. The number of possible attractors depends 
on the complexity of this dynamics and on the number 
of genes involved. It is widely believed that morpho- 
genesis is a precise, well-controled series of bifurcations 
which happen in the proliferating and migrating popula- 
tion of cells, generating an ever increasing complexity of 
patterns of differentiated cell regions. 

Disregarding some early asymmetric clevages, and non- 
uniform distribution of citoplasmic factors in the fertil- 
ized egg, cell divisions usually produce equivalent daugh- 
ter cells. Consequently, cell proliferation leads to an in- 
creasing domain of identical cells, where all the system 
parameters are distributed uniformly. Such a subsystem 
of identical cells is, however, embedded in, and com- 
municates with other, eventually already differentiated, 
groups of cells. 

There are essentially two ways that a subsystem of 
identical cells can later differentiate, either as a whole, 
or in parts: (i) There is a critical number of cells above 
which the spatially homogeneous attractor loses stabil- 
ity (Turing instability), leading to spontaneous spatial 
patterning. Q It is the size (the number of cells) of the 
increasing domain that plays the role of a bifurcation pa- 
rameter, (ii) There is an external inductive signal, eme- 
nating from an other group of already differentiated cells, 
which acts as a bifurcation parameter, and drives the sys- 
tem into the new, spatially inhomogeneous state. In this 



latter case, timing of the signal is crutial. 

Many biological examples could be mentioned for 
the two mechanisms of differentiation. Size-driven in- 
stabilities [case (i)] take place, e.g., in early insect 
development. A well-studied case is the syncytial blas- 
toderm stage of the fruit fly Drosophila, where a series of 
patterns of gene expression arise, forming various stripes 
of high and low concentration regions of gene products 
along the anterior-posterior axis. 

Inductive differentiation [case (ii)], on the other hand, 
is typical in later stages of development. Q As examples, 
we can mention the mesoderm and notochord induction 
in vertebrates, |l| the vulva formation in the soil nema- 
tode Caenorhahditis elegans, 5] and the development of 
the retina of the Drosophila fly,i^ where inductive influ- 
ence of the environment cells were clearly demonstrated. 

Our aim in this paper is to study a simple example of 
inductive differentiation. Emphasis will be put on the as- 
pects of cellular discreteness. The fact that interacting 
cells are discrete objects is usually overlooked in mod- 
elling biological pattern formation processes. However, 
as will be demonstarted here, spatial discreteness is a 
source of a variety of new phenomena with possible bio- 
logical significance. 



II. THE MODEL 

In the following we consider a semi-infinite one- 
dimensional chain of cells where the cell distance (lat- 
tice constant) is set to unity. We suppose that each cell 
in this system is characterized by the concentration of 
a single chemical (the morphogene) whose value (low or 
high) informs us about the actual state (type) of the cell. 
The morphogene concentration in cell n at time t will 
be denoted by Un{t). In an obviously highly oversimpli- 
fied setup the complicated cell biochemistry is reduced 
to an effective nonlinear autochatalitic reaction involv- 
ing the morphogene. We also assume that the morho- 
gene is diffusive and that the differentiation process can 
be described on a reaction-diffusion basis. Since the cells 
are discrete objects their diffusive coupling is modelled 
by a discrete Laplacian, and as it will be demonstrated. 
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this has far-reaching consequences. The analysis in the 
foUowing can be readily generalized to two- or three- 
dimensional domains with a straight surface if fluctua- 
tions in u{n, t parallel with the surface can be neglected. 

Inside the bulk of the system 2 < n < oo, our 
reaction-diffusion equation for the concentration distri- 
bution Un{t) takes the form 

= F{Un) + D{Un+l + Un-l - 2Un), (1) 

where D is the diffusion constant, and F{u) is a nonlinear 
reaction kinetics function characterising the cells in the 
bulk. We assume that the cell system is not coupled dif- 
fusively to its environment, but by receptor molecules in 
the cell membrane, it is capable of receiving an external 
inductive signal. Since real biological signal transduction 
mechanisms are complicated cascades of different cnzimc 
reactions, without worrying about the details here, we 
only assume that due to the signal the reaction kinet- 
ics function in the cells change. We consider the case 
when the penetration depth of the signaling molecules 
is so short that the signal is received almost exclusively 
by the very first cell along the line, and, for simplicity, 
the signal is thought to effect the morphogene production 
linearly. (Note that we make a clear distinction between 
the signaling molecules and the morphogene. The latter 
can freely diffuse in the system, while the former can- 
not.) With this proviso, we write the reaction-diffusion 
equation for the first cell in the form 

^ = S + Fim)+Diu2-ui). (2) 
ot 

Out of the various theoretical possibilities, in the fol- 
lowing we analize the case when F{u) is bistable and 
piece-wise linear 

, . / -f3u if w < a , , 

with (3 > and < a < 1. This is a charicature of the 
widely used Nagumo reaction kinetics function 

F{u) = Pu(l - u){u- a). (4) 

In the sequal /3 will be set /3 = 1, which can always be 
achieved by rescaling appropriately t, D and S. 

In both cases of F{u) the reaction kinetics of the cell 
is bistable in the lack of the external signal. When the 
signal is present, it acts as a chatalizer in cell 1 and in- 
creases the production rate of the morphogene. Even 
though the continuous form Eq. |0J is more realistic, we 
study in detail the piecewise linear charicature since it is 
analiticaly more tracktable. Numerical simulations car- 
ried out using the Nagumo form Eq. show that the 
qualitative behavior of the two models are essentially the 
same. Some minor differences will be pointed out in the 
sequel. 

In order to be able to assess the role of discreteness in 
the model we will also consider its usual continuous space 



analog, i.e., when the discrete Laplacian is replaced by 
the second derivative 

n ~> X 

du „d^u 

^-^^+^(-)- ->0 (5) 

Coupling to the environment via the S term in Eqs. (O 
translates into a Neumann boundary condition at the 
surface of the system du{x,t)/dx\x=o ~ S. 

The discrete and continuum models only become 
equivalent in the large D limit. This can be easily 
shown by dimensional analysis. The only parameter 
whose dimension contains the spatial length is D, [D] — 
[rn^/s]. As [(3] = [l/s] any solution u{x,t) of the con- 
tinuum model must contain x and D in the combination 
x/y/D/p. When D is large u{x,t) varies slowly in space 
so the second derivative can be discretized on the lattice 
without committing much error. Note that the discrete 
version contains an additional length scale: it is the lat- 
tice constant which was chosen to be unity. The solution 
of the discrete model is expected to deviate considerably 
from that of the continuous model when the diffusion 
length becomes comparable to the lattice constant, i.e.. 

The set of equations defined above contains the basic 
elements to model cellular differentiation in response to 
an external signal: Before switching on the inductive sig- 
nal, our system is uniform (undifferentiated). The mor- 
phogene concentration in every cell is u„ = 0, which is 
clearly a stable steady state. We can say that all the cells 
have type 0. When the external signal begins to increase 
(that we suppose to be adiabatically slow), the first cell 
at the end of the chain goes through a bifurcation, and 
switches from the branch u < a (type-0) to the branch 
u > a (type-1). It becomes differentiated. As the signal 
strength increases further, more and more cells flip into 
type 1. This avalanche of bifurcations may become self- 
sustaining, and the differentiation may sweep through the 
system in the form of a travelling wave. Under different 
conditions, the position of the domain wall separating 
type-0 and type-1 cells stays a well-defined function of 
the signal strength S. Then a natural question is what 
happens when S (adiabatically) returns to its original 
zero value. (According to biological observations, induc- 
tive signals are only present in a certain time interval of 
the process of development.) As we will see soon, eventu- 
ally the already differentiated cells do not de-differentiate 
into type 0, but maintain their type-1 state even in the 
lack of external signal. The system shows hysteresis. 

III. PROPAGATION FAILURE 

We begin our analysis with the classification of the 
possible bulk (i.e., far from the surface) behaviors. We 
analyze under what conditions can a two-domain steady 
state solution exist, when u„ is a monotonic decaying 
function of the cell position n and u^oo = 1, i*oo = 0. We 
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suppose that the domain wall (kink) is located between 
sites M — 1 and M, so that 

> a if n < M - 1 , , 

<a if n > M, 

where we introduced the superscript M to explicitely de- 
note the position of the kink. A concentration distribu- 
tion Un{t) satisfying Eq. © will be called a kink-M. 

It is well-known 7] that in the continuum version of 
the model in Eq. ^ for an infinite system (S* = 0), a 
steady state kink can only exist in the special case when 
a = 1/2. If a > 1/2, a kink- type initial profile develops 
instead into a travelling wave in which the domain wall 
travels with a constant speed c leftward. On the other 
hand, if a < 1/2, the doman wall travels rightward. 

In the lattice version of Eq. however, steady state 
domain wall solutions persist in a wide range of a values. 
When a e u^], with u± = u±{D), the domain wall is 
pinned and its propagation is impeeded. This, so called, 
propagation failMre|8ll9L llClj| is due to spatial discreteness. 
Travelling wave behavior only exits when a < U- or a > 
u+. 

In the case of the piece- wise linear function F{u) in Eq. 
the calculation of u- and u+ is straightforward. ;9j A 
candidate kink-M steady state solution of Eq. ^ with 
dun/dt = can be looked for in the form 

M _{ 1 + Ae"« if n < M 

Substituting this ansatz into Eq. , the inverse diffusion 
length At and the two constants A and B turn out to be 

K = 2sinh"^ V'TT^D (8) 



This allows us to identify the pinning region boundaries 
for a given D as 

«± = i±i(l + 4D)-V2. (12) 

The values of u+ and u_ are plotted in Fig. ^ as a func- 
tion of the diffusion constant D. Clearly, the obtained 
kink-M steady state solution is only valid in the shaded 
region of the diagram. On the other hand, when (D, a) 
is outside the shaded domain, the solution in Eqs. Q- 
is only a spurious solution. In this region of (£>, a) 
there are no steady state solutions; all initial conditions 
develop into travelling waves. 

Even though the effect of lattice pinning is alike for 
the case of the Nagumo-type (continuous) reaction func- 
tion, exact calculation of the pinning boundaries is not 
feasible. A perturbative app roach in the small a limit 
was carried out in Ref. [iflj. There is also an impor- 
tant difference how the wavefront speed c scales as, for 
a given a, the diffusion canstant approaches its critical 
value Dc = Dc{a). Simple bifurcation theory analvsisfioj 
shows that in the continuous F{u) case c scales following 
a power law with an exponent 1/2, i.e., 

c^iD^D.y/^ (13) 

while for the piece-wise linear kinetics the singularity is 
logarithmic 

-l/\n{D- Dc). (14) 

The latter form arises essentially from the nonanalicity 
(jump discontinuity) of the F(u) function in Eq. |(2Jl, and 
has been analysed in detail in Ref. 9]. 

IV. INDUCTION WITH HYSTERESIS 
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Note however that when the ansatz in Eq. iQ is used, 
one tacitly assumes that all cells on the left (right) of 
the kink are on the high (low) concentration branch of 
the piece-wise linear function F{u). Having found the 
solution in Eqs. (|8I9|I . this assumption must be checked 
for consistency: The concentration values obtained for 
the left (L = M - 1) and right (i? = M) neighboring 
cells of the kink arc 



,M 



= - + - tanh - = - + -(1 + 4D) 
22 2 22^ ' 



-1/2 



< = i-itanh| = l-l(l + 4Z?)-V2, (10) 

thus the above calculation is only consistent if we find 
that 



> a and 



u¥ < a. 



(11) 



Let us know investigate the inductive situation. As the 
signal S increases the distribution of concentration values 
Un{t) in the system becomes monotonically decreasing 
as a function of n. Even when the system is not in an 
equilibrium state we can define an M value characterizing 
the actual position of the domain wall separating type-1 
and type-0 cells using Eq. When seeking a steady 
state kink at site M, we solve the semi- infinite set of 
equations defined by Eqs. and Q with dun{t)/dt — 0. 
Again, the analytic solution is attainable for the piece- 
wise linearized kinetics. Working with the Ansatz 
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the unknown coefficients turn out to be 
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with K again given by Eq. Unlike in the translation 
invariant (infinite chain) case in Eq. (|1U|I . the concentra- 
tion values at the kink, ufj and , are now explicite 
functions of the kink position M and the signal strength 
S 

u'JiS)=e-^u^'. (19) 
In the limit M ^ oo we get back the bulk results in Eq. 

As it was done for the infinite system, consistency of 
the solution must be checked at this point. When the 
consistency condition of Eq. fails no steady state 
solution with the kink at site M exists. As a consequence 
the process of differentiation cannot stop at site M, and 
the domain wall moves on. 

Since the explicit expression for u'^ {S) and (S) is 
available in Eq. H19|) . for any values of /3, D, a and S we 
can readily construct the set of possible AI values {M} 
for which the consistency condition in Eq. holds, 
and thus the kink-M steady state exists. Although in 
theory every element of this set {M} could be realized 
as a steady state, it is the previous history of the system 
(the initial conditions) and the dynamics of the reaction- 
diffusion process which determines which steady state (if 
any) gets finally realized. This is in contrast with the 
continuum space model description where the position of 
a steady-state kink is always uniquely determined by the 
actual model parameters. 

The set of possible steady states on the a vs plane 
for the case /3 = 1,D — 2 fixed is depicted in Fig. [21 The 
different domains are separated by straight lines which 
is an artifact stemming from the simple form of F{u) in 
Eq. (PJ. Nevertheless, a qualitatively similar diagram can 
be obtained using the continuous Nagumo form. There 
are three main possibilities for a given a and 5*: (i) the 
number of steady states is finite, (ii) any M yields a valid 
steady state above a certain value (this is indicated by a 
"-I-"), (iii) there are no steady state kinks at all. 

In the piece-wise linear model under investigation the 
kink steady states, if they exist, are always stable against 
perturbations. Thus if at a given time t the actual kink 
position is not an element of the set of steady states {M}, 
differentiation or de-differentiation continues untill the 
domain wall reaches the first M value which is already 
in {M}. Having reached the domain of attraction of a 
stable steady state the kink stabilizes at that point and 
the process halts untill, eventually, a further change in S 
destabilizes the system again. 

Let us consider now an adiabatically slow process in 
which the inductive signal increases from zero to S'max 
and then decreases back to zero again. Using the above 
rule we can easily construct the phase diagram shown in 
Fig. for such a process. Domains are labelled by the 
M values of the kinks as they get realized in order. For 
example, the small domain [01230] has a history in which 
M increases continuously from to 3, then as the signal 



diminishes it jumps ubruptly back to 0. Once again there 
are three main regions, separated by thick lines, in this 
phase diagram: 

(i) In the upper part of the phase diagram the sys- 
tem gradually differentiates and then completely de- 
differentiates as the signal varies. The de-differentiation 
process can be continuous or may contain sudden jumps 
when the value of a is closer to u+. Note that this re- 
gion corresponds more or less to the values of a where 
in the infinite model kinks develop into travelling waves 
moving leftwards, i.e. towards the surface of our semi- 
infinite system. Due to this bias a continuous presence 
of the signal is needed to maintain differentiated type-1 
cells in the system. 

(ii) In the middle part of the phase diagram, correspond- 
ing approximatively to the pinning region of Fig. ^ the 
cells remain differentiated even when S falls back to zero. 
Note that when a is close to u_ the domain wall can take 
a huge jump at the beginning as S reaches a certain value. 
In this situation the maximum value of the signal S'max is 
important, since this is the factor which determines the 
range of the irreversibly differentiated domain. 

(iii) Finally for small values of a the differentiation pro- 
cess becomes self-sustaining when a critical value of the 
signal is exceeded. This mimics the infinite chain be- 
haviour with a travelling wave moving rightwards. The 
signal only triggers the differentiation process but after 
that it plays no further role. 

Typical examples of the three kinds of behavior are de- 
picted in Fig. ^a-c). 

The structure of the phase diagram in Fig. 13 is rather 
involved, demonstrating that even a simple model like 
this can show an amazing complexity. When the more 
realistic Nagumo-type reaction function is considered the 
exact solvability of the problem is lost. Nevertheless, 
numerical simulations we carried out demonstrated that 
the main conclusions about the qualitative behavior of 
the three phases remain unchanged. 



V. SUMMARY AND DISCUSSION 

In this paper we analysed a semi-infinite one- 
dimensional one-chemical reaction-diffusion system. The 
reaction kinetics was assumed to be bistable, giving rise 
to two different type of cells: type-0 (low concentration 
type) and type-1 (high concentration type). Starting 
from a homogeneous situation (all cell are type-0) the 
system underwent a differentiation process in response 
to an external inductive signal. The signal was modelled 
as a boundary condition in the continuum version and as 
an extra term in the internal cell kinetics of the first cell 
in the discrete space version. 

Depending on the model parameters the differentiation 
process sweeps through the whole system, or flips a lim- 
ited number of cells to type-1 up to a given position. We 
found that the behavior of the system differs consider- 
ably in the continuum and in the discrete space versions. 
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In the former the position of the domain wall between 
the two cell types is either a well-defined function of the 
external signal strength, or the front of differentiation 
inevitably develops into a travelling wave. In the dis- 
crete case the fate of the system depends on its previous 
history, giving rise to hysteresis. 

We analysed in detail the situation when the bistable 
reaction function is piece-wise linear. The model was 
solved analytically, and we constructed a detailed phase 
diagram based on the different types of behavior as a 
function of the model parameters. We found three mayor 
scenarios for the system: (i) In response to the inductive 
signal the solution develops into a travelling wave which 
differentiates the whole (semi-infinite) domain, (ii) the 
signal causes some spatially limited differentiation but 
when it diminishes all cells de-differentiate, (iii) differen- 
tiated cells get stabilized and the inhomogenious solution 
persists even when the signal disappears. 

Although the analysis was carried out with a some- 
what special reaction function, numerical simulations we 
have done support our expectation that the observed be- 
havior is widely universal in discrete space models, and 
the qualitative chategories found remain valid in simi- 
lar models with more realistic reaction functions. There 
are, of course, minor quantitative differences such as the 
type of scaling near the bifurcation points, or the actual 
domain wall location for a given signaling scheme. 

In general we have found that adding spatial discrete- 
ness to reaction diffusion models has a tendency to im- 
prove domain wall stability between different tissues, and 
to make the emerging pattern less susceptible to fluctua- 
tions of the signaling mechanisms. Since robostness and 
stability of developmental processes and that of the adult 
organism is an inevitable necessity for the survival of bi- 
ological species, we may wonder that the invention of 
cellular membranes by evolution, which made the fun- 
damental building blocks discrete, was at least in parts 
motivated by such a developmental benefit, among oth- 



ers. 

Finally we would like to mention an actual biological 
observation which seems to be explainable on the basis 
of the above model. During the retina differentiation 
of Drosophila it has been observed that ommatidia (the 
basic functional units of the retina consisted of photore- 
ceptor and other types of cells) develop behind a slowly 
moving wave front, the " morphogenetic furrow" . (6(1 There 
are many genes which are only expressed behind the fur- 
row, and thus one (or more) of them is believed to play 
the role of a morphogene. It was also noticed^J that a 
slight shift in the environmental temperature is enough 
to make the wave front stop. Untill the temperature is 
raised back to normal again the process of differentiation 
does not continue. This slight artificial manipulation, al- 
though capable of empeeding the propagation of the front 
for hours or days, is believed to have no residual effects 
on further retina development. 

Knowing that the propagating wavefront is extremely 
slow, c «?m/s, we can speculate that the developing 
retina is in fact tuned very close to a pinning region 
boundary. A change in the tissue temperature neces- 
sarily alters the actual model parameters. Although it 
would be very difficult to estimate on a phenomenolog- 
ical basis how the complex, non-linear set of biochemi- 
cal reactions get modified by a temperature decrese, we 
can at least assume that the overall diffusivity of the 
chemicals get reduced. This can drive the system into 
the pinning region as is illustrated in Fig. figrbulk, caus- 
ing eventually a propagation failure. In this unfavorable 
temperature regime the domain wall, represented by the 
morphogenetic furrow, between undifferentiated and al- 
ready differentiated cells becomes a stable steady state, 
and the differentiation process temporarily halts. 
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Fig. 1 
Path & Domanski 



FIG. 1: Phase diagram on the a vs D plane for an infinite chain. Pinned steady state solutions exist in the shaded region, 
while traveling waves exist in the unshaded ones. The pinning transition takes place along the a = a+{D) and the a = a-{D) 
curves, and can be initiated either by changing D for a fixed, or by varying a with D fixed (see arrows). 
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Fig. 2 
Path & Domanski 



FIG. 2: Steady state kink solutions on the a vs S plane with D = 2. The set {M1M2 ■ ■ ■} denotes the possible positions of the 
kinks in the given region. A "+" represents that all M values are possible above the preceding value. In the shaded regimes 
many small phases appear. 
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Fig. 3 
Path & Domanski 



FIG. 3: Phcise diagram for the signaling scheme 5 : — > Smax — > on the a vs Smax plane with D = 2. The symbol [M1M2 • • •] 
denotes the order of kink positions as they get realized. In the shaded regime many small phases appear. 
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Fig. 4 
Path & Domanski 



FIG. 4: Schematic hysteresis diagrams for the signaUng scheme 5:0^ Smax — > showing the actual domain wall position 
M as a function of the signal strength S. (a) all cells de-differentiate, (b) some cells remain differentiated, (c) all cells become 
differentiated as a travelling wave emerges. 



